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I Abstract. We describe a statistical model to estimate the covariance matrix of matter tracer 

^ two-point correlation functions with cosmological simulations. Assuming a fixed number of 

^^ cosmological simulation runs, we describe how to build a 'statistical emulator' of the two- 

00 

r — point function covariance over a specified range of input cosmological parameters. Because 

t^^ the simulation runs with different cosmological models help to constrain the form of the 

■^ covariance, we predict that the cosmology-dependent covariance may be estimated with a 

^^ comparable number of simulations as would be needed to estimate the covariance for fixed 

cosmology. Our framework is a necessary first step in planning a simulations campaign for 

analyzing the next generation of cosmological surveys. 
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1 Introduction 

Cosmological large-scale structure statistics contain valuable information about many cos- 
mological parameters. Accurate parameter inference from large-scale structure observations 
requires a model for the sample variance distribution of the observed statistics, for example 
the uncertainties and correlations between angular scales in the galaxy two-point correlation 
function. Often, the sample variance distribution is assumed to be multivariate Gaussian, 
requiring specification of only a covariance matrix of the cosmological two-point functions. 
While the covariance of the two-point function is known for a Gaussian field, effects of survey 
masks, galaxy clustering bias, and nonlinear gravitational evolution require many simulated 
realizations of the survey via A-body codes to accurately predict the covariance structure. 
Ref. [1] recently showed that errors in covariance estimates propagate into increased cos- 
mological parameter errors, which have the same effect as a reduction in the survey area. 
Near-term surveys will require at least W^ simulation realizations to prevent effective survey 
area losses larger than 10%. But, the computational requirements for covariance estimation 
will be even more challenging than forecasted in Ref. [1] because the cosmology dependence 
of the covariance is likely to become important as cosmological parameter constraints shrink 
in future surveys [2, 3]. 

In some analyses, the covariance matrix was estimated from the data using resampling 
methods such as the Jackknife or bootstrap. This has two large disadvantages. First, [4] 
showed that a variety of resampling methods underestimate both the variances and correla- 
tions compared to those derived from many simulation realizations. Second, any data-derived 



covariance estimator necessarily ignores the cosmology-dependence in the covariance which 
can significantly bias results [2]. 

In this paper we aim to determine how many cosmological simulations are required to 
achieve a target uncertainty in the covariance of the two-point correlations of mass density 
tracers. We also show how the number of simulation runs can be reduced over the brute force 
approach by exploiting the smooth variation of the covariance components as a function of 
standard cosmological parameters in a simulation emulator. 

We briefly argue for the need to model the cosmology dependence of the covariance in 
Section 2, although this has also been established in the literature [2, 3]. In Section 3 we 
specify the statistical model for the cosmology-dependent covariances that allows estimation 
of covariances with a minimal number of simulations. In Section 4 we apply the Fisher 
matrix to forecast the uncertainties in covariance matrix elements as functions of the number 
of cosmological simulations run. We summarize our main conclusion in Section 5 that it 
is more efficient to simultaneously model the CDC than run many simulation realizations 
at several fixed cosmological models. We also briefly describe the halo model used in this 
analysis and assess its accuracy the Appendix. 

2 Why cosmology-dependent covariances? 

To motivate a two-point statistic covariance emulator, it is helpful to understand how cos- 
mological parameter inference is affected by the cosmology dependence of the covariance. 
First, there is an issue of data interpretation: When comparing a model to data is it more 
helpful to ascribe uncertainties to the data (i.e. from observational limitations) or the model 
(i.e. from the statistical formulation of the statistics that are observed, commonly referred 
to as 'sample variance')? Second, we should determine how biased our parameter constraints 
could be if we erroneously ignore the cosmology dependence of the covariance. This latter 
issue was already addressed in [2] who showed that the inferred parameter constraints can 
change by many standard deviations in some cosmic shear scenarios and can also be biased. 
We leave similar predictions using our covariance emulator to future work. 

As stated previously, estimating the covariance model from the data tends to underes- 
timate the uncertainties [4] when compared with ensembles of simulated observations. Simu- 
lating the error on the model is then the preferred method, however, the model error is itself 
often model dependent. 

For example, consider A^ observed samples Xi from a one-dimensional zero-mean Gaus- 
sian distribution where we aim to estimate the standard deviation of the distribution. The 
maximum likelihood estimator for the standard deviation is. 
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This estimator has an error. 
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var(a^) = -:^^, (2.2) 

so the error on the estimator is also dependent on the value of the model parameter o"^ . This 
is a simple example, but the same principle holds for many cosmological estimators. 

We show an example of parameter-dependent errors on the angular correlation function 
in Figure 1 when the parameter is the square-root of the amplitude of the two-point function. 
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Figure 1. A cartoon example illustrates how parameter-dependent model errors can be important. 
A given set of observations, orange points, may be well fit by a model with a small amplitude and 
smaller model errors (green shaded region) or by a larger amplitude model that also has larger model 
errors (purple shaded region). If the errors from the model with the smaller amplitude are erroneously 
assumed to hold for all models we get the magenta band around the higher amplitude model; leading 
to a worse model fit to the orange data points. 

(Tg. In this example, the data points are generated with erg = 2 (for illustration only) and are 
consistent with models that have erg = 1.85 or 2.25 when the proper error is used. However, 
if one were to erroneously assign the errors for ug = 1.85 to models with varying cjg in 
the mean correlation function (i.e. a model-independent error), the model with the larger 
true variance would be incorrectly excluded by the data. This situation always arises when 
assuming fixed errors for cosmological models with different cg (although at a less obvious 
level for observationally consistent values of ug < 1). The same principle applies to inference 
of other cosmological parameters as well. 



3 Covariance matrix emulator and simulation design 

The covariance matrices of large-scale structure probes, C are typically estimated by run- 
ning many A^-body simulations with different pseudo-random number seeds in the initial 
conditions and constructing a sample covariance estimator from the outputs. This is a com- 
putationally intensive task, requiring 10^ A^-body simulations to reduce the errors in the 
covariance below other systematic uncertainties for current surveys and several fold for sur- 
veys of the future. These simulations must also be run for different cosmologies in order to 
properly give unbiased constraints. We address both of these issues with a statistical emula- 
tor described in this section that allows for a reduction in the number of simulations required 
as well as linking together the simulations run at different cosmologies. 
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simulations run at na points in parameter space with Ur^i 



independent simulation realizations at each point. For each simulation a summary statistic 
y^ (such as the power spectrum or correlation function) is computed (where k = 1, . . . , Ur^i). 



Given the set of y^, statistics we now describe how to estimate a model for the cosmology- 
dependent covariance (CDC) by means of a 'simulation emulator'. 

Simulation emulators have been successfully developed to model the mean matter power 
spectrum over a six-dimensional cosmological parameter space at high precision [5] and have 
been applied successfully to problems where the input parameters are less well-constrained [6- 
9]. 

An emulator is specified in 2 steps. The 'simulation design' defines at which points in 
cosmological parameter space A^-body simulations will be run to calibrate the emulator. The 
'emulation' step consists of calibrating a statistical model to interpolate the outputs of the 
simulation design runs to new points in parameter space. 

The first step in building a simulation design is choosing the points in parameter space 
where simulations will be run. The Orthogonal Array Latin Hypercube has proven to be a 
successful algorithm for choosing design points (but see also the improvements in [10]). 

Given specified design points, we then have to decide how many simulation realizations 
to run at each design point. We could require many simulation realizations at each design 
point if we need to construct a converged sample covariance estimator. However, with a 
careful parameterization of the covariance, we can use the simulation runs at all design 
points simultaneously to jointly constrain the covariance model at each design point. 

We achieve our goal in several steps with four features of the covariance emulator that 
all contribute to improve the estimate of the CDC over a brute-force sample covariance 
estimator, 

1. a careful decomposition of the two-point covariance matrix, 

2. a specification of orthogonal basis functions to decorrelate the covariance's components, 

3. optimized Gaussian Process parameters derived from the ensemble of simulation runs, 

4. calibration of the emulator to constrain the decorrelated mode amplitudes of the co- 
variance matrix components. 

Briefly, we first decompose our matrices using the Generalized Cholesky Decomposition to 
produce values that are easily emulated. Second, we create basis functions from these values 
using Principal Component Analysis (PGA). Third we fit (via maximum- likelihood) Gaus- 
sian Process parameters to link together the different cosmologies in the simulation design. 
Finally, we constrain our PGA mode amplitudes using the Gaussian Process. 

An important consideration that we neglect in this paper is whether a (potentially noisy) 
sample covariance estimator is first needed at each design point before the covariance emulator 
can be constructed. We avoid this issue by using an analytic model for the covariances, but 
Ref. [11] showed that sample covariance estimates might not be needed at any stage in the 
calculation. We will explore this further in a later publication. 

3.1 Covariance matrix decomposition 

We use the Generalized Cholesky Decomposition (GCD) to decompose the covariance ma- 
trices [following 12-14], which is alternately written in the following forms, 

C = LDL^ (3.1) 

C-^ = T^D-^T, T = L-\ (3.2) 



where T is a lower-triangular matrix with ones on the diagonal and D is a diagonal semi- 
positive definite matrix. The primary utility of the GCD for our purposes is that a positive 
definite covariance matrix, C, is guaranteed as long as all diagonal entries of D are positive, 
for any real values (j) in the lower triangular components of T. We can therefore interpolate 
unconstrained values of (p and In D over the cosmological parameter design space. 

We will index the GCD components at each simulation design point with i = 1, . . . , rid- 
Assuming the two-point function has n^ bins, we label the non-trivial components of T as [12] , 

T;., = -(/.;, ior2<j<nb,£ = l,...,j-l. (3.3) 

Similarly, we label the diagonal components of D as, 

D;.^. = exp(4) fori <j<n6. (3.4) 

3.2 Basis functions 

The GCD provides a method to separate the nfc(nf,+l)/2 unique components of the covariance 
matrix into a computationally convenient form. We can further improve the covariance 
estimation by reducing the number of components that must be modeled as functions of 
cosmology. At the same time, we can use all the simulation runs, with all cosmological 
inputs, to determine common structures in the covariance matrices. 

Following the emulator construction in [5, 11], we proceed by extracting the components 
of InD and T into separate vectors for each design point, stack these vectors for all design 
points into an Uy x 71^ array (where Uy is either equal to n^ or nh{nh — l)/2), and compute 
Principal Components (PC) to identify a subset of components that are most strongly varying 
over the design space. We choose to compute principle components of In D so that we obtain 
values supported on the entire real line. 

Following [10] we first subtract the mean of the design run components d and (p and 
then scale the result by a single number so that the combined entries of our Uy x n^ matrix 
of design run components have variance one. 

d^ = add' + dc (3.5) 

ct>' = cj^4>' + 4>c- (3.6) 

We collectively label the centered and scaled matrices by /x = d or 0, and perform a singular 
value decomposition: /x = UBV where U has dimension ny x p {p = min(ny,nd)) with 
U"^U = Ip, V has dimension na x p with V'^V = Ip, VV^ = I„^, and B {p x p) is a diagonal 
matrix of singular values. The matrix of basis vectors, = -t^^UB with weights w = y^nrfV 



/rid 

normalized so that —ww = I„^ (this latter choice makes it simple to specify priors on the 
Gaussian Process variance parameters). 

We keep only the first Pd,P<i> ^ P columns of 0. Then we redefine the scaled design 
covariance matrix components as a sum over pD or p,^ modes, 

pd 

i=i 

P-P , 

4^tm = Y,5}^^,tj, t = l,...,-nbinb-l) (3.7) 



where ^x,tj is the entry of the matrix <Px in row t (t = 1, . . . , n^) and column j {j = 1, . . . ,p£, 
or p^), ■^j{0) and 5j{6) are the jth (parameter-dependent) mode ampHtudes (y^V-^ above), 
and 6 are the model (i.e. cosmological) parameters. Ref. [11] included i.i.d. Normal errors in 
the truncated mode decompositions to account for residuals when PD,(f) < P (which will always 
be assumed so that the PCA achieves some reduction in parameters). For the purposes of 
forecasting we ignore the error in the truncation of the PCA expansion, but this error should 
be propagated when constructing a full emulator of the covariance (see e.g. Eqs. 7-11 of 

[11])- 

Note that [12] suggest a similar decomposition for the simultaneous modeling of several 

covariance matrices. Our Eq. 3.7 differs in using basis vectors that are independent of the 

model parameters 6 while imparting all the model dependence to the mode amplitudes in 

the truncated basis. 

Joint estimation of the covariance matrices is now reduced to estimating the mode 
amplitudes 7^ and (5*- from the simulations run at each design point, given the basis vectors 
0£) and 00, which are estimated from the combination of all simulation runs. 

We define estimators for the mode amplitudes from the likelihood for the simulation 
design runs given the model in Eq. 3.7. Ref. [12] show that with the GCD, the log-likelihood 
for the set of covariance matrices at the simulation design points can be written, 



rid nt, 



^(y^|7J,^i) = EE ir'f^ - ^^^A.exp(-d^)) , (3. 



=1 i=i 
where Tjj is the jth column of Tj and 

' k=l 

is the sample covariance matrix estimate at design point i. The mean y* could either be 
specified by a theoretical model (e.g. halofit) or could be the sample mean from the 
simulation realizations at each design point. Ref. [11] consider how the sample mean can 
be jointly estimated with the sample covariance, but here we assume the mean has zero 
uncertainty. 

We can then estimate the mode amplitudes conditioned on the simulation design runs 
either with maximum-likelihood estimators or posterior samples. In order to use simulation 
runs with different input cosmologies to jointly constrain a covariance model, we now must 
specify how the mode amplitudes 7^(0) and 5j{9) can vary with 6. 

3.3 Parameters for the Gaussian Process 

We impose Gaussian Process (GP) priors on the mode amplitudes 7*- and (5* as a means of 
linking the covariance models with different input cosmologies (in addition to the common 
structure imposed by the basis vectors •$>£) and 0,^) as well as a means of interpolating 
the covariance predictions to regions of parameter space where no simulations have been 
run [5, 10, 11, 15]. The GP priors impose requirements on the smoothness of the mode 
amplitude surfaces over the cosmological parameter space. 

We model each mode amplitude as independent GPs with zero mean (because we already 
subtracted the mean values of the d and (f) over the simulation design) and an exponential 
covariance model described by a single precision parameter for each mode Xxj and correlation 



parameters for each direction in the pg-dimensional cosmofogical parameter space px,ij £ 

[0,1] [11], 

ne,e',px,,\x,) = ^,\{pS;'^ , (3.10) 

where X = D,(f) and j = 1, . . . ,Pd,(I>- The correlation parameters control the smoothness of 
the mode amplitudes along each parameter direction. If all px,ij ~ 1 then the modes of the 
covariance are highly correlated over the design parameter space, allowing all simulations to 
aid in the estimation of the joint covariance model. Conversely, when many of the correlation 
parameters are near zero, we expect the covariance emulator to give little improvement 
beyond the standard sample covariance estimators at each design point. 

Our argument for using the covariance model in Eq. 3.10 rather than some other model 
is simply because it works, as we show in Section 4. The Matern covariance is a more flexible 
covariance model that is often used for GPs. Our choice of the covariance model in Eq. 3.10 
is partly informed by our expectation that the cosmological covariances we are interested 
in will be smoothly varying over the parameter spaces that are already tightly constrained 
by existing CMB and large-scale structure observations. It is possible that parameters in 
alternative cosmological models (e.g. modified gravity or dynamical dark energy) will be less 
constrained and a more flexible GP covariance model could be appropriate. 

Restricted to the design points 6*, the priors on the design mode amplitudes are then 
multivariate Normal, 

Trhj\Xd,j,Pd,j) = iV(0, S(r ; prf,„ A^,,)) (3.11) 

7r(d,|A^j,p^j) = N{0,^{e*;p^j,X^j)), (3.12) 

and the posterior for forecasting constraints on the mode amplitudes given the design runs 
is obtained by combining Eqs. 3.8 and 3.11, 

Phj^ '^ily*> ^d,j, Pd,j, X^,j, P<t>,j) = L {yl\Yj, S'j)'^{lj\>'d,j, pd,j)TT{Sj\X^j,p^,j). (3.13) 

We later compute the Fisher matrix for the mode amplitudes by taking derivatives of Eq. 3.13. 

The statistical model for the CDC is completed by calibrating the parameters of the 
GP models for each mode "fj and Sj. So, the simulation design runs are used to both infer 
the structure of the covariance parameterization in the form of the basis vectors <i>D and 
0(^ and to infer the cosmology dependence in the form of the GP parameters for the mode 
amplitudes. 

Previous emulators in the cosmology literature built hierarchical Bayesian models to 
marginalize over the GP parameters for each mode amplitude and thereby propagate all 
interpolation and parameterization uncertainties. We have built a similar framework here 
(in Eq. 3.13), but for the purposes of forecasting, will now use maximum-likelihood estimators 
for the GP parameters. 

By definition, the likelihood for the Gaussian process model for a mode amplitude 
evaluated at the design points is multivariate Gaussian [eq. (5.8) of 16], 

ln(p(y|a)) = -^y^S'iy - ^ In |S,| - ^ln2vr, (3.14) 

where y = <Sj;i = l,...,nd\ or < 7*; i = 1, . . . , n^^ > as determined by our fiducial model, 
and we have jointly labeled the GP model parameters as a = {p. A}. Taking derivatives with 



respect to the GP parameters, [16] find (their eq. 5.9), 

— ln(p(y|a)) = -Trf (aa^-S-i)— j witha = S-V (3.15) 

We also include hyperpriors on the parameters of the GP models, 

Px 

vr(Ax) = n^xy'e-^^^^- (3.16) 

i=i 

Px P0 

So we use the same hyperprior parameters for all modes of d and all modes of (f). 

3.4 Mode amplitude constraints 

Next we use the Fisher matrix to forecast the errors on the components of the covariance 
matrix estimates at every point in our simulation design given the set of simulation design 
runs. We aim to determine how many cosmological simulations are needed to achieve a given 
precision in the elements of the covariance matrix estimates given: 

• the number of simulation design points n^, 

• the number of simulation realizations at each design point n^^j, 

A related question is whether it is better to estimate the high-precision covariances at 
a few fixed cosmological models, or to estimate the CDC jointly with a few realizations at 
many cosmological parameter values. The answer to this question will partly depend on 
how smoothly the chosen components of the covariance matrices vary over the cosmological 
parameter space. We will show that the GCD and PCA parameterization of the model 
covariances described in section 3 yield such smoothly varying parameters for a standard 
wCDM model. 

The Fisher matrix is defined as the negative of the expectation of the curvature of the 
log-likelihood (or log-posterior as given in Eq. 3.13) about its peak, 

'dH{y,e_ 

^'j I e=eo 
The terms contributing to the Fisher matrix for the CDC mode amplitudes are 
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d'^D,k^D,£ 



2 l2 
= 9^r,i— ^(^M (3-19) 

2 rid 

where $fc is the kth row of the covariate matrix ^, bj is the jth singular value in the PC 
decomposition of the design runs, and we have used the definition of <& as composed from the 
SVD of the design runs. So, the conditional Fisher information for the 7 mode amplitudes 
is the same for all design points, is independent for each mode amplitude, and only depends 



on the SVD of the design runs (in the form of the variance of each mode over the design). 
Eq. 3.19 is analogous to the standard error on the variance. 
For the 5 modes however, 



d^ 



d6i^d6j^ 



n. 



i^lTr (tri ($<^,^) C tri ($^,^)^ D-^) , (3.20) 



where $0^,^ is the column of <t>^ indexed by rj and tri (<&</,,,;) indicates we force ^^^rj to fill 
the lower-triangular elements of an nf, x rif, matrix with all other entries zero. So the Fisher 
information depends on the covariance O at each design point i. And, the covariance between 
different 6 modes is nonzero, which means that constraining some modes of (p (at a fixed design 
point) can help constrain the other modes as well. 

The joint Fisher matrix for 7 and 6 has nonzero cross-terms, 






d5id-i\ 



tri ($<^,^) (r) ^ diag ($D,m) + diag ($z),m) (T') ^ tri ($<^,^) 



(3.21) 

Notice that eqs. 3.19, 3.20, and 3.21 scale linearly with n^^j as might be expected. From 

this we can immediately see that without a model for connecting the mode amplitudes at 

different design points, it is not possible to reduce the error in the covariance matrix elements 

faster than x/ru-- 



4 Simulation design study: cosmic shear 

To asses the performance of the covariance matrix emulator, we use the analytical halo 
model [17] to predict covariance matrices of the shear correlation function. Previous studies 
have shown that the halo model captures the correct qualitative features of the nonlinear 
two-point function covariances. While our modeling is expected to lack precision relative to 
what would be estimated from A^-body simulations, we believe the complexity of the model is 
sufficient to demonstrate the utility of our statistical framework. Also, any plans for running 
a large number of simulations to estimate the CDC must initially rely on imperfect models. 
For details on the model used see Appendix A. 

We assume the same cosmological parameters and similar ranges of variation as in Ref . [5] , 
shown in Table 1. All examples assume 32 design points (i.e. n^ = 32) in an Orthogonal 
Array Latin Hypercube (OALH) spanning this 5- dimensional design space. At each of the 
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0.611 


1.011 
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0.86 


1.06 
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-0.7 



Table 1. Ranges of the cosmological parameters for our example OALH simulation design. 

32 design points we compute a model for the nonlinear covariance of the shear correlation 
function assuming a 5-function source distribution at z = 1 and negligible shape noise. 
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Figure 2. The range of values of the covariance of the shear correlation function over the simulation 
design space. Left: diagonal terms in the covariance. Right: rows of the matrix of correlation 
coefficients. Each panel represents a different bin in angular scale in the correlation function (in 
arcmin.). For the right pannel, we see the qualitative trends we expect, strong correlation on small 
scales and weaker at large scales. The bands with larger widths at arcminute scales is due in part to 
the change in angular diameter distances with cosmology. 
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Figure 3. Distribution of components of the Generalized Cholesky Decomposition of the design 
covariance matrices for the different cosmology design points. In this decomposition, a majority of 
the variation in the matrix elements is confined to a few indices. 
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Figure 4. The first 8 principal components of tlie simulation design ln(_D) (black circles) and cj) (red 
triangles) components of the Generalized Cholesky Decomposition of the shear correlation function 
covariance matrix. We keep all PC modes that contribute a fraction of more than 10^"^ to the total 
variance (i.e. all modes above the horizontal dashed line). 

4.1 Covariance matrix decomposition 

In Figure 2 we show the full range of values in covariance diagonal components (left) and 
correlation coefficients (right) over the design space. The diagonal terms of the covariance 
span well over an order of magnitude in value and differ in shape as the angular projection of 
the one- halo term changes. The amount of cross-correlation is mostly stable with cosmology 
however, it does vary greatly at large scales. Figure 3 shows the amount of variation in the 
GCD compoments, d and cf) over the full simulated space. In this decomposition we see much 
of the variation in the covariance has been confined to fewer components. 

4.2 Basis functions 

In Figure 4 we show the first 8 PCA amplitudes as a function of mode index. For both the 
diagonal elements d and lower triangular elements (f) we see that most of the variation is 
contained within the first few modes of the decomposition allowing for significant reduction 
in the dimensionality of the covariance matrices over the full simulation design. We retain 
those PC modes that contribute at least 10~^ to the fractional variance. 

4.3 Parameters of the Gaussian Process 

As stated previously in Section 3.4 the amount of correlation in the GP modes informs us 
how strongly dependent the information in the CDC is on different cosmological parameters. 
Figure 5 shows the mode amplitude correlations (i.e. the px parameters in Eq. 3.10) with 
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Figure 5. Maximum likeliiiood estimates for the GP correlation parameters for the first 3 PC modes 
for the diagonal (left) and lower-triangular (right) components of the GCD of the covariance. A 
correlation parameter close to one indicates the mode amplitude is smoothly varying along a given 
parameter axis. Conversely, a correlation parameter near zero indicates large variations in the mode 
amplitude along the given axis. 

each cosmological parameter in a 'sensitivity' analysis. A value of one means that the mode 
amplitudes are highly correlated along a given parameter axis, so the emulator is not sensitive 
to variations in this parameter. This also means sparser sampling of simulation design runs 
can be used. The parameters wq, ilbh'^, and Ug are often the most correlated, indicating they 
have little impact on the covariance. The mode amplitudes are consistently vi^eakly correlated 
along the ag and O.mh'^ parameter axes, indicating these parameters largely determine the 
form of the covariance as might be expected because the amplitude of the shear correlation 



function depends on the product 



~ (t; 






4.4 Mode amplitudes of the covariance matrices 

The marginal Fisher errors of the mode amplitudes for the diaginal and lower triagular 
compoments of the GCD, 7 and 6 respectively, are shown in Figure 6 as functions of n^ 
for n^ = 32. Here we compare the emulator error performance in the mode amplitudes to 
that of a brute force simulation. The errors shown are for a single cosmology design point, 
however the emulator is constrained over the full parameter space considered. We see that 
the emulator performs as well or better depending on the number of simulation run and with 
the added benefit of modeling the CDC over the full parameter space. While the fractional 
error at higher modes is constrained by the GP prior, keep in mind that the prior is informed 
from the data. The GP prior can be thought of as the maximum variation over all the 
cosmology design points. 

If we compare the two lines in Fig. 6 for fixed fractional error (i.e. reading horizontally), 
we see that in many cases the number of simulations required for the emulator is about one 
order of magnitude smaller. The results for other n^ values are similar, except that the GP 
prior is less significant for a fixed rir- 
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Figure 6. Forecasted errors on the mode amplitudes of the diagonal (left) and lower-triangular (right) 
components in the GCD of the covariance matrices. The panels show the first 3 principle components 
in each case. The blue bands show the range of marginal errors over 32 design points when using all 
design points simultaneously to constrain the mode amplitudes. The orange bands show the range of 
marginal errors when the simulations run with different cosmologies are used The errors then scale 
with Y^rv- The dashed horizontal line indicates the variance in the GP prior imposed on each mode 
amplitude (where the GP prior parameters are determined by maximum-likelihood estimates given 
the model covariance matrices at each design point). 



We have not propagated any uncertainties from the calibration of the GP parameters 
or the truncation of the PC basis expansion. Including these uncertainties will tend to bring 
the emulator forecasts closer to those assuming no relationship between modes at different 
design points. But, given the success of other emulator frameworks in the literature, we do 
not expect additional sources of uncertainty to qualitatively change the results in Fig. 6. 



5 Conclusions 

Estimating covariance matrices of power-spectra and correlation functions is a computation- 
ally intensive task that requires many CPU hours of simulations to achieve the accuracy 
required for future surveys [18, 19]. In this paper we have shown that the precision in the es- 
timated cosmology dependent covariance matrices can only improve as fast as 1/y^, where 
rir is the number of simulation realizations, if one considers each cosmological model indepen- 
dently. If however, one simultaneously models the simulations at different cosmologies using 
Gaussian Processes, the number of simulations required to reach a given precision can be 
reduced while also modeling the CDC (which cannot otherwise be done with a sample of dis- 
joint cosmological simulations). This makes the computational challenges of simulating the 
analysis much more tractable. However, we do not find orders of magnitude improvements 
in the number of simulations needed and other methods of estimating the covariances should 
be considered in combination with the emulator presented here (e.g. shrinkage estimators 
[20], large-scale mode-resampling [21], and optimized simulation design spaces [10].) 

Future work will demonstrate the accuracy of the full emulator framework, with all 
uncertainties propagated including those from the lack of sample covariance estimators at 
each design point. 
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A The halo model 

For our 'simulation' of covariances we employ the halo model [17] to estimate the nonlinear 
matter power spectrum to be used in the analysis. This modeling is part of the python 
cosmology prediction package CHOMP^. We use the halo model as defined in [22] using a 
[23] mass function. This model is accurate to within 20% at 1 arcminute when compared to 
A^-body simulations from [24] and within 10% between 20-200 arcminutes. For the covari- 
ances we follow the formalism of several papers [25-27] , considering two contributions to our 
covariance matrix, a Gaussian and a non-Gaussian term. 

c{ei,e2) = CG{eue2) + CNGieue,) (a.i) 

For the Gaussian term we use the definitions as laid out in [26] using Limber's approx- 
imation 



cf\eue2) = ^jidiJo{iei)uie^){p^\i)p^\i) + p'\i)pi\i) + ^{5,1,5,1 + 5a6,k)] 

(A.2) 
where A is the area of the survey, a"^ is the variance per galaxy pair, n* is the density of 
galaxies for probe z, and P^^{1) is the projected power spectrum written as 

P''{1) = / dx ^^^^V -P(- ■ X) (A.3) 

Jo X ^ X 

where P is the non- linear matter power spectrum as a function of redshift, x is the comving 
distance, /*(x) is the weighted window function (e.g. the lensing kernel), and the integration 
is from to the horizon. 

For the non-Guassian term, we model only the one halo term of the halo-trispectrum. 
Following [28], this term is 

Ti-halo{h,k2)= [dM^{^)Uy{h,Mf*y{k2,Mf (A.4) 

J dM p 

Where ^ iwthe number density of halos as a fucntion of mass M, p is the matter 
density, and y is the Fourier transform of halo density profile normalized to 1 at k = 0. For 
this analysis we assume that halos follow an NFW [29] profile. Here we only considering 
the tri-spectrum in a parallelogram configuration as required by the covariance estimation. 
Projecting this configuration with the window functions we have 

V'iluh) = / r{x)f{x)fHx)f{x)/x-''T{^, - ■■ X) (A.5) 

Jo XX 

The final form of the non-Gaussian covaraince is 

C*Jfc'(0i,^2) = ^^ j hdh j l2dl2Jo{hOi)Jo{he2)r'''\hM) (A.6) 

It is worth noting that more terms contribute to the covariance and are highly dependent 
on the specified survey geometry. One parameterization as presented by [30] gives the beat 
coupling and halo sample variance as a single Super-Sample Covariance term. This term is 

^available at: http://code/google.coin/p/choinp 
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Figure 7. Comparison of CHOMP lensing two-point functions with those from [31], and [24] using 
ray-tracing through iV-body simulations. Left: Comparison between the simulations and the chomp 
shear correlations. Right: Comparison bewteen the simiulated and predicted covariance matrices. 
The correlation function is predicted to within 10% at acminute scales in the halo model and to 
sub-arcminute scales for HaloFit. The covariances agree to within 20% for scales larger than 10 
arcminutes. The results of this anlaysis hold as long as the variation of both the simulated and 
predicted covariances are smooth as a function of cosmology. 

as dominant as the non-Gaussian term at scales A; > 1 Mpc/h, however, we do not consider 
them in this analaysis currently. While this term is not sub-dominant, within the scope of this 
paper, we do not expect it to significantly change the covariance's dependence on cosmology 
or add any challenge to the marix decomposition and emulation presented in this work. 

B Validation of halo model shear correlation covariance 

To apply our forecasts for planning simulation runs, it is important to understand how well 
our covariance model may reproduce the covariance derived from simulations. In this section 
we compare our covariance model from the CHOMP code to the 1000 lensing simulation re- 
alizations from [24, 31]. These simulations give the power spectra and real-space correlations 
of cosmic shear for a variety of source redshifts. We compare our model constructed from 
the CHOMP cosmology package to these sims. The left pannel of Figure 7 shows the level of 
agreement between our models and the Sato results. We see good agreement between these 
models and the A^-body sims. For the halofit model derived in [32] we find agreement to 
within 5% at arcminute scales for the real-space correlation. The Halo Model used in this 
analysis fare's slighly worse but is still within 10% for similar scales to the HaloFit model. 

The right pannel of Figure 7 shows how well the halo model recovers the features of 
the Sato covarianes. Overall, we find agreement between the simluations at various scales. 
The reader should keep in mind, though, that the results of this paper hold as long as the 
covariance matrix variation with respect to cosmology is smooth and therefore still gives an 
adiquate description of the covariance as a funciton of cosmology. Future work will explore 
the robustness and limitations of this code further. 
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